home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / Scientist's Helper src / s.helper.5 / regseries.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-04-19  |  15.3 KB  |  631 lines  |  [TEXT/EDIT]
  1. #include "all.h"
  2. #include "regtabext.h"
  3.  
  4. histogram()
  5. {
  6.     int i, n, iCol, *x, *y;
  7.     float vMin, vMax, s, *u, binNumber, nm1;
  8.             
  9.     SToR( command.cmdWord[1], &vMin, FALSE);
  10.     SToR( command.cmdWord[2], &vMax, FALSE);
  11.     if( vMax<=vMin ) {
  12.         ErrMsg("bounds out of order");
  13.         }
  14.     SToI( command.cmdWord[3], &n );
  15.     if( (n<1) || (n>table.header.maxRows) ) {
  16.         ErrMsg("bad number of bins");
  17.         }
  18.     s=(vMax-vMin)/n;
  19.     nm1 = (float)(n-1);
  20.         
  21.     /* grab temporary memory */
  22.     x = NewPtr( (long)(2*n) );
  23.     if( MemError()!=noErr ) {
  24.         ErrMsg("not enough memory");
  25.         }
  26.  
  27.     for( iCol=2; iCol<=table.header.cols; iCol++ ) { /*loop over columns in table*/
  28.     
  29.         /*lock position of column */
  30.         HLock( table.ptr[iCol-1] );
  31.         
  32.         /* zero counters */
  33.         y=x;
  34.         for( i=1; i<=n; i++ ) {
  35.             *y=0;
  36.             y++;
  37.             }
  38.         
  39.         /* count up hits */
  40.         u=*(table.ptr[iCol-1]);
  41.         for( i=1; i<=table.header.rows; i++ ) {
  42.             binNumber = (*u-vMin)/s;
  43.             if( (binNumber>=0.0) && (binNumber<=nm1) ) {
  44.                 *(x+(int)binNumber)+=1;
  45.                 }
  46.             u++;
  47.             }
  48.             
  49.         /* copy to table */
  50.         u=*(table.ptr[iCol-1]);
  51.         y=x;
  52.         for( i=1; i<=n; i++ ) {
  53.             *u = (float)(*y);
  54.             u++;
  55.             y++;
  56.             }
  57.  
  58.         /*unlock current column*/
  59.         HUnlock( table.ptr[iCol-1] );
  60.         } /*next column*/
  61.         
  62.     /*dispose of temporary memory*/
  63.     DisposPtr( x );
  64.  
  65.     table.header.start = vMin;
  66.     table.header.samp = s;
  67.     table.header.rows = n;
  68.     table.header.interpolated=TRUE;
  69.     CreateCol1();
  70.     Header2Vars();
  71. }
  72.  
  73. taper()
  74. {
  75.     int xCol, outputCol, startRow, endRow, i;
  76.     float f, x, x1, x2, y, y2;
  77.     
  78.     if( (strcmp(command.cmdWord[1],"ascending")!=0) && (strcmp(command.cmdWord[1],"descending")!=0) ) {
  79.         ErrMsg( badModifier );
  80.         }
  81.     SToI( command.cmdWord[2], &xCol );
  82.     if( GoodCol(xCol)!=0 ) {
  83.         ErrMsg("bad x column");
  84.         }
  85.     SToI( command.cmdWord[3], &outputCol );
  86.     if( (GoodCol(outputCol)!=0) || (table.header.interpolated&&(outputCol==1)) ) {
  87.         ErrMsg("bad output column");
  88.         }
  89.     SToI( command.cmdWord[4], &startRow );
  90.     if( (startRow<1) || (startRow>table.header.maxRows) ) {
  91.         ErrMsg("bad starting row");
  92.         }
  93.     SToI( command.cmdWord[5], &endRow );
  94.     if( (endRow<1) || (endRow>table.header.maxRows) ) {
  95.         ErrMsg("bad ending row");
  96.         }
  97.     if( endRow<=startRow ) {
  98.         ErrMsg("start row, end row out of order");
  99.         }
  100.     GetTable( startRow, xCol, &x1 );
  101.     GetTable( endRow, xCol, &x2 );
  102.     f = 3.1415926 / (x2-x1);
  103.     
  104.     if( strcmp(command.cmdWord[1],"ascending")==0 ) {
  105.         for( i=startRow; i<=endRow; i++ ) {
  106.             GetTable(i, xCol, &x );
  107.             GetTable(i, outputCol, &y );
  108.             y2 = y * 0.5 * (cos( (double)(3.1415926+f*(x-x1)) ) + 1.0 );
  109.             SetTable(i, outputCol, y2, FALSE );
  110.             } /*end for i*/
  111.         } /*end if ascending*/
  112.     else if( strcmp(command.cmdWord[1],"descending")==0 ) {
  113.         for( i=startRow; i<=endRow; i++ ) {
  114.             GetTable(i, xCol, &x );
  115.             GetTable(i, outputCol, &y );
  116.             y2 = y * 0.5 * (cos( (double)(f*(x-x1)) ) + 1.0 );
  117.             SetTable(i, outputCol, y2, FALSE );
  118.             } /*end for i*/
  119.         } /*end if descending*/
  120. }
  121.  
  122. convolve()
  123. {
  124.     int inputCol, opCol, n, outputCol, i, j, k, nm1;
  125.     float *x, *u, *v, *w, sum;
  126.     
  127.     if( !table.header.interpolated ) {
  128.         ErrMsg("table must be interpolated");
  129.         }
  130.     SToI( command.cmdWord[1], &inputCol);
  131.     if( GoodCol(inputCol)!=0 ) {
  132.         ErrMsg("bad input column");
  133.         }
  134.     SToI( command.cmdWord[2], &opCol);
  135.     if( GoodCol(opCol)!=0 ) {
  136.         ErrMsg("bad operator column");
  137.         }
  138.     SToI( command.cmdWord[3], &n );
  139.     if( (n<1) || (n>table.header.rows) ) {
  140.         ErrMsg("bad operator length");
  141.         }
  142.     SToI( command.cmdWord[4], &outputCol);
  143.     if( (GoodCol(outputCol)!=0) || (outputCol==1) ) {
  144.         ErrMsg("bad output column");
  145.         }
  146.     
  147.     /* grab temporary memory */
  148.     x = NewPtr( (long)(4*table.header.rows) );
  149.     if( MemError()!=noErr ) {
  150.         ErrMsg("not enough memory");
  151.         }
  152.         
  153.     /* lock position of columns */
  154.     HLock( table.ptr[inputCol-1] );
  155.     HLock( table.ptr[opCol-1] );
  156.     HLock( table.ptr[outputCol-1] );
  157.     
  158.     /*do convolution*/
  159.     u=x;
  160.     v=*(table.ptr[inputCol-1]);
  161.     w=*(table.ptr[opCol-1]);
  162.     nm1=table.header.rows-1;
  163.     for( i=0; i<table.header.rows; i++ ) {
  164.         sum = 0.0;
  165.         for( j=0; j<n; j++ ) {
  166.             k = i-j;
  167.             if( (k>=0) && (k<=nm1) ) {
  168.                 sum += (*(v+k)) * (*(w+j));
  169.                 } /*end if*/
  170.             } /*end for j*/
  171.         *u=sum;
  172.         u++;
  173.         } /*end for i*/
  174.     
  175.     /*copy results to table*/
  176.     u=x;
  177.     v=*(table.ptr[outputCol-1]);
  178.     for( i=0; i<table.header.rows; i++ ) {
  179.         *v=(*u) * table.header.samp;
  180.         u++;
  181.         v++;
  182.         } /*end for i*/
  183.         
  184.     /* lock position of columns */
  185.     HUnlock( table.ptr[inputCol-1] );
  186.     HUnlock( table.ptr[opCol-1] );
  187.     HUnlock( table.ptr[outputCol-1] );
  188.  
  189.     /*dispose of temporary memory*/
  190.     DisposPtr( x );
  191. }
  192.  
  193. spect()
  194. {
  195.     int ifft, n, iCol, i;
  196.     float *x, *u, *v, *w;
  197.     
  198.     if( (strcmp(command.cmdWord[1],"amplitude")!=0) &&
  199.         (strcmp(command.cmdWord[1],"power")!=0) &&
  200.         (strcmp(command.cmdWord[1],"phase")!=0) ) {
  201.         ErrMsg( badModifier );
  202.         }
  203.     if( !table.header.interpolated ) {
  204.         ErrMsg("table must be interpolated");
  205.         }
  206.     
  207.     n = table.header.rows;
  208.     if( n<=128 ) { /*find nearest power of two*/
  209.         ifft=128;
  210.         }
  211.     else if( n<=256 ) {
  212.         ifft=256;
  213.         }
  214.     else if( n<=512 ) {
  215.         ifft=512;
  216.         }
  217.     else if( n<=1024 ) {
  218.         ifft=1024;
  219.         }
  220.     else if( n<=2048 ) {
  221.         ifft=2048;
  222.         }
  223.     else {
  224.         ifft=4096;
  225.         }
  226.         
  227.     /* grab temporary memory */
  228.     x = NewPtr( (long)(4*(ifft+2)) );
  229.     if( MemError()!=noErr ) {
  230.         ErrMsg("not enough memory");
  231.         }
  232.  
  233.     for( iCol=2; iCol<=table.header.cols; iCol++ ) { /*loop over columns in table*/
  234.     
  235.         /* lock position of column */
  236.         HLock( table.ptr[iCol-1] );
  237.         
  238.         /* swap into temporary array */
  239.         u=x;
  240.         v=*(table.ptr[iCol-1]);
  241.         for( i=1; i<=table.header.rows; i++ ) {
  242.             *u = *v;
  243.             u++;
  244.             v++;
  245.             }
  246.         for( i=table.header.rows+1; i<=ifft; i++ ) { /*zero rest of array*/
  247.             *u = 0.0;
  248.             u++;
  249.             v++;
  250.             }
  251.         
  252.         cfftr( x, ifft ); /*fast fourier transform*/
  253.         u=x;
  254.         for( i=1; i<=ifft+2; i++ ) { /*proper normalization*/
  255.             *u = (*u) * table.header.samp;
  256.             u++;
  257.             }
  258.         
  259.         if( strcmp(command.cmdWord[1],"amplitude")==0 ) {
  260.             u=x;
  261.             v=x+1;
  262.             w=*(table.ptr[iCol-1]);
  263.             for( i=1; i<=ifft/2+1; i++ ) {
  264.                 *w = (float) sqrt( (double)((*u)*(*u)+(*v)*(*v)) );
  265.                 u=u+2;
  266.                 v=v+2;
  267.                 w++;
  268.                 }
  269.             }
  270.         else if( strcmp(command.cmdWord[1],"power")==0 ) {
  271.             u=x;
  272.             v=x+1;
  273.             w=*(table.ptr[iCol-1]);
  274.             for( i=1; i<=ifft/2+1; i++ ) {
  275.                 *w = (*u)*(*u)+(*v)*(*v);
  276.                 u=u+2;
  277.                 v=v+2;
  278.                 w++;
  279.                 }
  280.             }
  281.         else if( strcmp(command.cmdWord[1],"phase")==0 ) {
  282.             u=x;
  283.             v=x+1;
  284.             w=*(table.ptr[iCol-1]);
  285.             for( i=1; i<=ifft/2+1; i++ ) {
  286.                 *w = (float) atan2( (double)(*v), (double)(*u) );
  287.                 u=u+2;
  288.                 v=v+2;
  289.                 w++;
  290.                 }
  291.             }
  292.  
  293.         /*unlock current column*/
  294.         HUnlock( table.ptr[iCol-1] );
  295.         } /*next column*/
  296.         
  297.     /*dispose of temporary memory*/
  298.     DisposPtr( x );
  299.  
  300.     table.header.start = 0.0;
  301.     table.header.samp = 1.0 / (ifft*table.header.samp);
  302.     table.header.rows = ifft/2+1;
  303.     CreateCol1();
  304.     Header2Vars();
  305. }
  306.  
  307. cfour( data, n, isign ) /*fast fourier transform on complex data*/
  308. int n, isign;
  309. float data[];
  310. {
  311.         int ip0, ip1, ip2, ip3, i3rev;
  312.         int i1, i2a, i2b, i3;
  313.         float sinth, wstpr, wstpi, wr, wi, tempr, tempi, theta;
  314.         ip0=2;
  315.         ip3=ip0*n;
  316.         i3rev=1;
  317.         for( i3=1; i3<=ip3; i3+=ip0 ) {
  318.                 if( i3 < i3rev ) {
  319.                         tempr = data[i3-1];
  320.                         tempi = data[i3];
  321.                         data[i3-1] = data[i3rev-1];
  322.                         data[i3] = data[i3rev];
  323.                         data[i3rev-1] = tempr;
  324.                         data[i3rev] = tempi;
  325.                         }
  326.                 ip1 = ip3 / 2;
  327.                 do {
  328.                         if( i3rev <= ip1 )
  329.                                 break;
  330.                         i3rev -= ip1;
  331.                         ip1 /= 2;
  332.                         }
  333.                         while ( ip1 >= ip0 );
  334.                 i3rev += ip1;
  335.                 }
  336.                 ip1 = ip0;
  337.                 while ( ip1 < ip3 ) {
  338.                         ip2 = ip1 * 2;
  339.                         theta = 6.283185/( (float) (isign*ip2/ip0) );
  340.                         sinth = (float) sin( (double) (theta/2.) );
  341.                         wstpr = -2.*sinth*sinth;
  342.                         wstpi = (float) sin( (double) theta );
  343.                         wr = 1.;
  344.                         wi = 0.;
  345.                         for ( i1=1; i1<=ip1; i1+=ip0 ) {
  346.                         for ( i3=i1; i3<ip3; i3+=ip2 ) {
  347.                                 i2a=i3;
  348.                                 i2b=i2a+ip1;
  349.                                 tempr = wr*data[i2b-1] - wi*data[i2b];
  350.                                 tempi = wr*data[i2b] + wi*data[i2b-1];
  351.                                 data[i2b-1] = data[i2a-1] - tempr;
  352.                                 data[i2b] = data[i2a] - tempi;
  353.                                 data[i2a-1] += tempr;
  354.                                 data[i2a] += tempi;
  355.                                 }
  356.                                 tempr = wr;
  357.                                 wr = wr*wstpr - wi*wstpi + wr;
  358.                                 wi = wi*wstpr + tempr*wstpi + wi;
  359.                                 }
  360.                         ip1=ip2;
  361.                         }
  362.                 return;
  363.                 }
  364.  
  365. cfftr( x, n ) /*fast fourier transform on real data*/
  366. int n;
  367. float x[];
  368. {
  369.         int nn, is, nm, j, i;
  370.         int k1j, k1i, k2j, k2i;
  371.         float s, fn, ex, wr, wi, wwr, wrr, wwi, a1, a2, b1, b2;
  372.         nn = n/2;
  373.         is = 1;
  374.         cfour( x, nn, is );
  375.         nm = nn/2;
  376.         s = x[0];
  377.         x[0] += x[1];
  378.         x[n] = s - x[1];
  379.         x[1] = 0.0 ;
  380.         x[n+1] = 0.0;
  381.         x[nn+1] = (-x[nn+1]);
  382.         fn = (float) n;
  383.         ex = 6.2831853 / fn;
  384.         j = nn;
  385.         wr = 1.0;
  386.         wi = 0.0;
  387.         wwr = (float) cos( (double) ex );
  388.         wwi = (float) (-sin( (double) ex ));
  389.         for (i=2; i<=nm; i++) {
  390.                 wrr = wr*wwr-wi*wwi;
  391.                 wi = wr*wwi+wi*wwr;
  392.                 wr = wrr;
  393.                 k1j = 2*j-1;
  394.                 k1i = 2*i-1;
  395.                 k2j = 2*j;
  396.                 k2i = 2*i;
  397.                 a1 = 0.5*(x[k1i-1]+x[k1j-1]);
  398.                 a2 = 0.5*(x[k2i-1]-x[k2j-1]);
  399.                 b1 = 0.5*(-x[k1i-1]+x[k1j-1]);
  400.                 b2 = 0.5*(-x[k2i-1]-x[k2j-1]);
  401.                 s = b1;
  402.                 b1 = b1*wr+b2*wi;
  403.                 b2 = b2*wr-s*wi;
  404.                 x[k1i-1] = a1-b2;
  405.                 x[k2i-1] = (-a2-b1);
  406.                 x[k1j-1] = a1+b2;
  407.                 x[k2j-1] =  a2-b1;
  408.                 j -= 1;
  409.         }
  410.         return;
  411. }
  412.  
  413. integrate()
  414. {
  415.     int i, xCol, yCol, outputCol;
  416.     float x1, x2, y1, y2, sum;
  417.             
  418.     SToI( command.cmdWord[1], &xCol);
  419.     if( GoodCol(xCol)!=0 ) {
  420.         ErrMsg("bad x column");
  421.         }
  422.     SToI( command.cmdWord[2], &yCol);
  423.     if( GoodCol(yCol)!=0 ) {
  424.         ErrMsg("bad y column");
  425.         }
  426.     SToI( command.cmdWord[3], &outputCol);
  427.     if( (GoodCol(outputCol)!=0) || ((outputCol==1)&&table.header.interpolated) ) {
  428.         ErrMsg("bad output column");
  429.         }
  430.         
  431.     sum = 0.0;
  432.     GetTable( 1, yCol, &y1 );
  433.     GetTable( 1, xCol, &x1 );
  434.     SetTable( 1, outputCol, sum, FALSE );
  435.     for( i=2; i<=table.header.rows; i++ ) {
  436.         GetTable( i, xCol, &x2 );
  437.         GetTable( i, yCol, &y2 );
  438.         sum = sum + 0.5*(y1+y2)*(x2-x1);
  439.         SetTable( i, outputCol, sum, FALSE );
  440.         x1=x2;
  441.         y1=y2;
  442.         }
  443. }
  444.  
  445. differentiate()
  446. {
  447.     int i, xCol, yCol, outputCol;
  448.     float sum, x1, x2, y1, y2;
  449.             
  450.     SToI( command.cmdWord[1], &xCol);
  451.     if( GoodCol(xCol)!=0 ) {
  452.         ErrMsg("bad x column");
  453.         }
  454.     SToI( command.cmdWord[2], &yCol);
  455.     if( GoodCol(yCol)!=0 ) {
  456.         ErrMsg("bad y column");
  457.         }
  458.     SToI( command.cmdWord[3], &outputCol);
  459.     if( (GoodCol(outputCol)!=0) || ((outputCol==1)&&table.header.interpolated) ) {
  460.         ErrMsg("bad output column");
  461.         }
  462.         
  463.     GetTable( 1, xCol, &x1 );
  464.     GetTable( 1, yCol, &y1 );
  465.     for( i=1; i<=table.header.rows-1; i++ ) {
  466.         GetTable( i+1, xCol, &x2 );
  467.         GetTable( i+1, yCol, &y2 );
  468.         SetTable( i, outputCol, (y2-y1)/(x2-x1), FALSE );
  469.         x1=x2;
  470.         y1=y2;
  471.         }
  472.     SetTable( table.header.rows, outputCol, infinity, FALSE );
  473. }
  474.  
  475. summation()
  476. {
  477.     int i, inputCol, outputCol;
  478.     float sum, x;
  479.             
  480.     SToI( command.cmdWord[1], &inputCol);
  481.     if( GoodCol(inputCol)!=0 ) {
  482.         ErrMsg("bad input column");
  483.         }
  484.     SToI( command.cmdWord[2], &outputCol);
  485.     if( (GoodCol(outputCol)!=0) || ((outputCol==1)&&table.header.interpolated) ) {
  486.         ErrMsg("bad output column");
  487.         }
  488.         
  489.     sum = 0.0;
  490.     for( i=1; i<=table.header.rows; i++ ) {
  491.         GetTable( i, inputCol, &x );
  492.         sum = sum + x;
  493.         SetTable( i, outputCol, sum, FALSE );
  494.         }
  495. }
  496.  
  497. bandpass()
  498. {
  499.     int i, inputCol, outputCol;
  500.     float flow, fhigh, digt, *fbout, *t, *v;
  501.     
  502.     if( !table.header.interpolated ) {
  503.         ErrMsg("cant bandpass uninterpolated table");
  504.         }
  505.         
  506.     digt = 1.0 / table.header.samp;
  507.     SToR( command.cmdWord[1], &flow, FALSE);
  508.     SToR( command.cmdWord[2], &fhigh, FALSE);
  509.     SToI( command.cmdWord[3], &inputCol);
  510.     if( (GoodCol(inputCol)!=0) || (inputCol==1) ) {
  511.         ErrMsg("bad input column");
  512.         }
  513.     SToI( command.cmdWord[4], &outputCol);
  514.     if( (GoodCol(outputCol)!=0) || (outputCol==1) ) {
  515.         ErrMsg("bad output column");
  516.         }
  517.     
  518.     /* grab temporary memory */
  519.     fbout = NewPtr( (long)(4*table.header.rows) );
  520.     if( MemError()!=noErr ) {
  521.         ErrMsg("not enough memory");
  522.         }
  523.      
  524.     /* bandpass input */
  525.     HLock( table.ptr[inputCol-1] );
  526.     t=*(table.ptr[inputCol-1]);
  527.     FILT( t, fbout, table.header.rows, 0, table.header.rows-1, digt, flow, fhigh, 2 );
  528.     HUnlock( table.ptr[inputCol-1] );
  529.     
  530.     /* swap result back into table */
  531.     HLock( table.ptr[outputCol-1] );
  532.     t=fbout;
  533.     v=*(table.ptr[outputCol-1]);
  534.     for( i=1; i<=table.header.rows; i++ ) {
  535.         *v = *t;
  536.         v++;
  537.         t++;
  538.         }
  539.     HUnlock( table.ptr[outputCol-1] );
  540.  
  541.     /*dispose of temporary memory*/
  542.     DisposPtr( fbout );
  543.     
  544. }
  545.  
  546.  
  547. FILT( FBIN, FBOUT, NN, K, M, DIGT, FLOW, FHIGH, N)
  548.  
  549.     float FBIN[], FBOUT[], FLOW, FHIGH, DIGT;
  550.     int NN, K, M;
  551.     {
  552.     
  553.     /* CHEBYSHEV RECURSIVE DIGITAL BANDPASS DIGITAL FILTER */
  554.     
  555.     /* FBIN:    INPUT ARRAY TO FILTERED */
  556.     /* FBOUT:    RESULTANT FILTERED ARRAY */
  557.     /* NN:        LENGTH OF FBIN, FBOUT */
  558.     /* K, M:    DATA FILTERED BETWEEN THESE LIMITS */
  559.     /* IDIGT:    SAMPLING RATE */
  560.     /* FLOW:    LOW PASS FREQ */
  561.     /* FHIGH:    HIGH PASS FREQ */
  562.     /* N:        ORDER OF FILTER (MUST BE EVEN OR IT WILL BE MADE EVEN) */
  563.     
  564.     float D[6], G[6], PI=3.1415926, ES=1.0, E, FL, FH, CF, BW, FAC, AYE, BEE;
  565.     float PIN, ARG, SR, SI, TEMP;
  566.     int N2, I, J;
  567.     
  568.     E=0.1; /* ten percent ripple */
  569.     if( (N&1)!=0 ) N=N+1; /*only even orders allowed*/
  570.     FL=2.0*FLOW/DIGT;
  571.     FH=2.0*FHIGH/DIGT;
  572.     
  573.     /* CALCULATE BANDWIDTH AND CENTER FREQUENCIES IN FRAME WARPED BY Z-FORM */
  574.     CF = 4.0 * tan( (double) (FL*PI/2.0) ) * tan( (double) (FH*PI/2.0) );
  575.     BW = 2.0 * ( tan( (double) (FH*PI/2.) ) - tan( (double) (FL*PI/2.) ) );
  576.     FAC = pow( sqrt((double)(1.0+1.0/(E*E))) + 1.0/E, (double)(1.0/N) );
  577.     AYE = 0.5*(FAC-1.0/FAC);
  578.     BEE = 0.5*(FAC+1.0/FAC);
  579.     PIN = PI/(2.0*N);
  580.     N2=N/2;
  581.     
  582.     for( I=1; I<=N2; I++ ) { /* LOOP THROUGH PAIRS OF CONJUGATE POLES */
  583.         ARG = (2*I-1)*PIN+PI/2.0;
  584.         SR = AYE * cos( (double)ARG);
  585.         SI = BEE * sin( (double)ARG);
  586.         ES = sqrt( (double)(SR*SR+SI*SI) );
  587.         TEMP= 16.0 - 16.0*BW*SR + 8.0*CF + 4.0*ES*ES*BW*BW - 4.0*BW*CF*SR + CF*CF;
  588.         D[1] = 1.0;
  589.         D[2] = 4.0*(-16.0 + 8.0*BW*SR - 2.0*BW*CF*SR + CF*CF)/TEMP;
  590.         D[3] = (96.0 - 16.0*CF - 8.0*ES*ES*BW*BW + 6.0*CF*CF)/TEMP;
  591.         D[4] = (-64.0 - 32.0*BW*SR + 8.0*BW*CF*SR + 4.0*CF*CF)/TEMP;
  592.         D[5] = (16.0 + 16.0*BW*SR + 8.0*CF + 4.0*ES*ES*BW*BW + 4.0*BW*CF*SR + CF*CF)/TEMP;
  593.         TEMP = 4.0*ES*ES*BW*BW/TEMP;
  594.         G[1] = TEMP;
  595.         G[2] = 0.0;
  596.         G[3] = -2.0*TEMP;
  597.         G[4] = 0.0;
  598.         G[5] = TEMP;
  599.         CONVL(FBIN,FBOUT,NN,K,M,D,G);
  600.         } /*end for i*/
  601. }
  602.  
  603.  
  604. CONVL(FBIN,FBOUT,NN,K,M,D,G)
  605.  
  606.     float FBIN[], FBOUT[], D[], G[];
  607.     int NN, K, M;
  608.     {
  609.     
  610.     int I, J, IX;
  611.     float TEMP;
  612.     
  613.     for( I=K; I<=M; I++ ) { /*SUM PAST VALUES OF INPUT*/
  614.         TEMP=0.0;
  615.         for( J=1; J<=5; J++ ) {
  616.             IX=I-J+1;
  617.             if( IX>=0 ) {
  618.                 TEMP=TEMP+FBIN[IX]*G[J];
  619.                 }
  620.             } /*end for J*/
  621.         
  622.         for( J=2; J<=5; J++ ) { /* SUM PAST VALUES OF OUTPUT */
  623.             IX=I-J+1;
  624.             if( IX>=0 ) {
  625.                 TEMP=TEMP-FBOUT[IX]*D[J];
  626.                 }
  627.             } /*end for J*/
  628.     
  629.         FBOUT[I]=TEMP;
  630.         } /*end for I*/
  631. }